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A  quasi-two-dimensional  numerical  model  is  presented  for  the  efficient  computation  of  the  steady- 
state  current  density,  species  concentration,  and  temperature  distributions  in  planar  solid  oxide  fuel  cell 
stacks.  The  model  reduction  techniques,  engineering  approximations,  and  numerical  procedures  used  to 
simulate  the  stack  physics  while  maintaining  adequate  computational  speed  are  discussed.  The  results 
of  the  model  for  benchmark  cases  with  and  without  on-cell  methane  reformation  are  presented  with 
comparisons  to  results  from  other  research  described  in  the  literature.  Simulations  results  for  a  multi¬ 
cell  stack  have  also  been  demonstrated  to  show  capability  of  the  model  on  simulating  cell  to  cell  variation. 
The  capabilities,  performance,  and  scalability  of  the  model  for  the  study  of  large  multi-cell  stacks  are  then 
demonstrated. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

A  fuel  cell  is  an  electrochemical  conversion  device  that  produces 
electricity  directly  from  oxidation  of  a  fuel.  Fuel  cells  are  charac¬ 
terized  by  their  electrolyte  material  and,  as  the  name  implies,  the 
solid  oxide  fuel  cell  (SOFC)  has  an  oxide  ceramic  electrolyte  [1]. 
Some  advantages  of  this  class  of  fuel  cell  include  high  efficiency, 
electrolyte  stability,  fuel  flexibility,  low  emissions,  high  power  den¬ 
sity,  simpler  fabrication  compared  to  other  cell  types,  and  the 
availability  of  byproduct  heat  for  cogeneration;  however,  the  great¬ 
est  disadvantage  is  the  high  operating  temperature  that  results 
in  mechanical/chemical  compatibility  issues,  thermal  stresses,  and 
longer  start-up  and  shut-down  times.  Fuel  cell  technology,  includ¬ 
ing  ceramic  electrolytes,  has  existed  for  well  over  a  century,  and 
industrial  research  into  fuel  cell-based  power  sources  has  been 
ongoing  since  the  late  1950s  [2]. 

Over  the  last  30  years,  research  teams  around  the  world  have 
shifted  more  focus  to  solid  oxide  fuel  cell  technologies  as  evidenced 
by  the  rapid  rise  of  research  publications  of  which  the  bulk  is 
related  to  SOFCs  [3].  In  the  United  States,  the  Solid  State  Energy 
Conversion  Alliance  (SECA)  was  formed  by  the  U.S.  Department 
of  Energy  in  2000,  and  partnerships  were  forged  among  industry, 
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national  laboratories,  and  academia  to  make  SOFCs  a  cost-effective 
alternative  power  source  [4-6].  As  part  of  this  expanding  research 
effort,  the  physical  phenomena  occurring  in  SOFCs  have  been  stud¬ 
ied  using  numerical  modeling  to  gain  scientific  understanding  and 
engineering  guidance.  Numerical  models  are  crucial  to  accelerating 
knowledge  and  understanding  within  any  SOFC  development  activ¬ 
ity  because  they  complement  the  expensive  and  time  consuming, 
but  necessary,  experimental  testing.  Fast  and  efficient  numerical 
modeling  can  be  used  to  perform  virtual  experiments  that  are 
important  to  cell  design  and  stack  operation,  thereby  reducing  the 
resources  and  development  time  needed  for  a  new  SOFC  concept. 

Models  for  studying  SOFCs  at  different  length  scales  continue  to 
be  important  to  advance  knowledge  and  understanding  of  how  to 
improve  performance,  stability,  reliability,  and  efficiency.  Molec¬ 
ular  and  microstructural  scale  models  have  been  important  for 
understanding  the  fundamental  electrochemical  behaviors  occur¬ 
ring  in  the  cell  and  the  effect  of  materials  and  structure  on  the  cell 
performance  (e.g.,  [7-1 1  ]).  Cell  and  stack  level  models  are  essential 
for  the  SOFC  designer  to  achieve  good  power  output  while  mini¬ 
mizing  temperatures  and  thermal  gradients  that  lead  to  structural 
degradation  or  failure.  Numerous  cell/stack  level  models  have  been 
generated  to  evaluate  the  coupled  flow,  electrochemical,  and  ther¬ 
mal  behaviors  of  various  stack  geometries  (e.g.,  [12-15,17-19]). 
SOFCs  are  desirable  for  power  generation  [16],  so  system-level 
models  are  used  to  obtain  the  operating  conditions  and  controls 
necessary  to  achieve  the  highest  possible  efficiency  for  the  SOFC  or 
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Fig.  1.  Illustration  of  (a)  co-flow/counter-flow  and  (b)  cross-flow  planar  SOFC. 


hybrid  generation  systems  based  on  the  unique  operating  require¬ 
ments  of  the  SOFC  module  [20].  More  extensive  reviews  of  the 
fundamental  mathematical  models  for  fuel  cells  and  developed 
SOFC  numerical  models  have  been  compiled  by  various  authors 
[16,21-25]. 

The  work  presented  in  this  paper  attempts  to  bridge  the  stack 
and  system  levels  by  bringing  internal  spatial  resolution  to  a 
computationally  feasible  submodel  for  use  in  system  level  stud¬ 
ies.  Other  researchers  have  also  recognized  the  need  to  move 
away  from  the  typical  zero-dimensional  (0D  or  ID,  2D,  or  3D 
for  one-,  two-,  or  three-dimensional)  thermodynamic  models  that 
provide  performance  based  solely  on  the  stack  average  state  to 
create  efficient  numerical  models  that  resolve  internal  distribu¬ 
tions  [26-28].  In  prior  work  by  the  authors,  tools  that  solve  the 
coupled  electrochemistry,  thermal,  and  flow  analysis  for  realistic 
3D  cell  geometries  were  created  in  both  the  computational  fluid 
dynamics  (CFDs)  and  finite-element  frameworks  [13,29].  While 
tubular  [30,31],  microtubular  [32-35],  flattened  tubular  [36],  and 
segmented-in-series  cells  [37]  are  being  pursued  by  multiple  teams 
and  remain  active  viable  concepts  for  SOFC  geometries,  the  authors 
have  focused  primarily  on  planar  designs,  which  exhibit  the  highest 
power  density. 

While  the  3D  model  is  extremely  useful  for  engineering  design 
of  the  cell  geometry  or  short-stack  performance,  the  computational 
requirements  for  detailed  modeling  of  modern  tall  stacks  increases 
dramatically.  Stack  developers  using  SOFCs  for  megawatt-scale 
power  generation  applications  are  striving  for  larger  and  taller 
stacks  to  reduce  costs  where  30-kW  modules  of  nearly  200  cells  in 
series  are  being  operated,  and  planar  areas  approaching  1000  cm2 
are  being  evaluated  [38,39].  To  achieve  the  high  efficiencies  for 
SOFC  or  hybrid  SOFC-turbine  power  plants  using  system  model¬ 
ing,  models  beyond  purely  thermodynamic  calculations  are  needed 
to  identify  peak  internal  temperatures  or  current  densities  impor¬ 
tant  to  stack  reliability.  Toward  this  end,  the  3D  numerical  model 
developed  previously  can  be  reasonably  simplified  to  a  more  com¬ 
putationally  efficient  2D  model  for  symmetric  flow  configurations. 
In  this  paper,  we  present  such  a  model  that  more  rapidly  computes 
the  solution  for  symmetric  multi-cell  stacks  while  still  providing 
spatial  resolution  for  the  temperature,  current,  and  species  distri¬ 
butions  along  the  direction  of  gas  flow. 

2.  Model  descriptions  and  assumptions 

Simulating  the  SOFC  stack  operation  requires  consideration  of 
multiple  physical  phenomena,  such  as  heat  transfer,  mass  transfer, 
charge  transfer,  chemical  reactions,  electrochemical  behavior,  and 
structural  deformation.  Since  fully  coupled  computational  models 
considering  all  of  the  above  fields  simultaneously  with  high  fidelity 


are  not  yet  available,  existing  SOFC  modeling  efforts  have  been 
focusing  on  specific  aspects  of  the  SOFC  operation  or  design  with 
available  computational  resources.  For  example,  prior  work  within 
the  authors’  modeling  group  developed  simulation  tools  to  evaluate 
the  electrochemical  performance  of  3D  cell  geometries  [13,29],  and 
the  predicted  thermal  fields  were  subsequently  used  for  mechani¬ 
cal  stress  analyses  to  assess  structural  performance  [39-41]. 

While  the  detailed  3D  model  will  provide  information  from  the 
stack-level  fuel/oxidant  outlet  temperatures  or  power  output  as  a 
function  of  stack  temperature,  its  computational  cost  is  usually  high 
for  solving  large  multi-cell  stack  towers  for  megawatt-scale  power 
generation  systems,  for  which  case  a  faster  modeling  tool  would  be 
more  desirable.  On  the  other  hand,  existing  system-level  models 
typically  offer  no  spatial  resolution,  resulting  in  a  0D  or  ID  result 
that  neglects  important  thermal  gradients  and  mischaracterizes  the 
temperature  state  of  the  stack. 

The  objective  of  this  effort  was  therefore  to  provide  an  efficient 
computational  tool  that  can  provide  the  users  with  distributions  of 
the  current  density,  voltage,  species  concentration,  and  tempera¬ 
ture  within  the  multi-cell  stack.  The  ability  of  the  user  to  define  all 
model  parameters,  material  properties,  boundary  conditions,  and 
current-voltage  relationship  is  still  required  as  well  as  the  capa¬ 
bility  to  study  the  effect  of  multiple  cells  stacked  in  series.  The 
developed  tool  has  proven  to  be  useful  in  characterizing  the  thermal 
state  of  SOFC  stacks  of  up  to  100  cells  with  computation  durations 
measured  in  minutes  rather  than  hours  or  days.  This  paper  dis¬ 
cusses  a  quasi-2D  SOFC  multi-physics  model  (SOFC-MP)  that  solves 
the  coupled  electrochemistry  analysis  and  heat  transfer  calcula¬ 
tions  in  a  2D  cross-section  of  the  stack  to  give  detailed  internal 
profiles  of  temperature,  current  density,  flow  compositions,  and 
cell  voltage  distribution.  This  is  based  on  the  user-defined  geome¬ 
try,  initial  flow  conditions,  SOFC  electrochemistry,  and  boundary 
conditions.  The  basis,  assumptions,  formulation,  and  numerical 
implementation  of  this  model  and  demonstrations  of  its  capabil¬ 
ities  are  described  in  the  following  sections. 

2.1.  Assumptions  for  the  quasi- two- dimensional  approximation 

Many  SOFC  manufacturers,  including  several  of  the  SECA  pro¬ 
gram  participants,  have  used  the  planar  configuration  [6],  so  this 
has  been  the  primary  geometry  of  interest  for  the  authors.  In  the 
traditional  planar  design,  the  fuel  and  air  channels  can  be  arranged 
in  co-flow,  counter-flow,  or  cross-flow  configurations  to  supply  the 
fuel  and  oxidant  in  the  most  efficient  manner  to  maximize  power 
and  minimize  temperature  gradients  within  the  stack  (Fig.  1 ).  A  3D 
model  is  required  for  detailed  evaluations  or  subsequent  structural 
analyses,  but  it  is  nevertheless  computationally  expensive.  FIow- 
ever,  when  only  co-flow  or  counter-flow  cases  are  considered,  the 
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Fig.  2.  Temperature  distribution  (°C)  in  the  active  cell  for  (a)  cross-flow,  (b)  co-flow,  and  (c)  counter-flow  configurations  [6]. 
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Fig.  3.  Current  density  distribution  (A cm-2)  in  the  active  cell  for  (a)  cross-flow,  (b)  co-flow,  and  (c)  counter-flow  configurations  [6]. 


3D  model  can  often  be  reasonably  approximated  by  a  ID  or  2D 
model  for  engineering  purposes  with  the  following  assumptions 
on  the  flow,  thermal,  and  electrochemical  behaviors. 

The  first  assumption  is  that  the  physical  fields  across  the 
width  direction,  perpendicular  to  the  flow  direction,  remain  largely 
unchanged.  For  example,  detailed  3D  results  of  Recknagle  et  al. 
[29]  presented  contour  plots  for  temperature,  species  concentra¬ 
tion,  and  current  density  distributions  over  the  cell  surface  with  the 
fuel  and  air  in  co-flow,  counter-flow,  and  cross-flow  configurations. 
All  of  the  computational  results,  including  the  temperature,  hydro¬ 
gen  concentration,  and  current  density  fields  shown  in  Figs.  2-4, 
are  quite  uniform  across  the  width  direction  for  the  co-flow  and 
counter-flow  cases.  Similar  observations  can  be  made  from  other 
modeling  results  [12,28,42-43].  Since  these  primary  results  of 
interest  show  little  variations  in  the  width  direction,  the  actual  3D 
problem  can  be  reasonably  approximated  by  a  2D  solution.  The 
assumption  becomes  less  appropriate  as  the  heat  loss  to  the  sides 
becomes  more  significant  or  as  the  cell  area  becomes  very  small. 


Nevertheless  the  assumption  should  be  generally  valid  for  large, 
well  insulated  cells  that  are  nearly  adiabatic. 

The  second  assumption  is  that  the  thickness  dimensions  of  a  cell 
are  small,  and  the  physical  properties  in  the  thickness  direction 
are  constant  for  each  domain  in  the  cell  (fuel  channel,  air  chan¬ 
nel,  interconnect,  and  the  positive  electrode/electrolyte/negative 
electrode  (PEN)  assembly).  For  the  solid  PEN  and  interconnect 
structures,  temperature  variation  through  the  thickness  is  assumed 
to  be  small.  For  fuel  and  air  gas  channels,  the  mass  transport  is 
treated  as  a  fully  developed  laminar  flow  described  sufficiently  by 
its  mean  velocity  and  temperature.  Therefore,  through-thickness 
gradients  in  these  domains  are  ignored,  but  temperature  differ¬ 
ences  between  these  domains  are  permitted  based  on  the  thermal 
transport  mechanisms.  This  assumption  leads  to  further  simplifi¬ 
cation  of  our  model  such  that  only  four  control  volumes  are  used 
along  the  out-of-plane  direction  of  the  cell,  e.g.,  interconnect,  fuel, 
PEN,  and  air.  Since  the  interconnect  serves  a  dual  purpose  as  the 
gas  barrier  and  the  current  collector,  it  often  has  a  complex  struc- 
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Fig.  4.  Hydrogen  mass  distribution  (%)  in  the  active  cell  for  (a)  cross-flow,  (b)  co-flow,  and  (c)  counter-flow  configurations  [6]. 
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Fig.  5.  Illustration  of  SOFC-MP  2D  mesh. 


ture  (e.g.,  ribs,  channels,  mesh,  stamping)  that  permits  flow  while 
transmitting  current  between  a  series  of  connected  cells.  In  SOFC- 
MP  2D,  effective  properties  are  used  instead  of  explicitly  modeling 
the  geometric  details  of  interconnect  to  include  their  influence  on 
conductive  heat  transfer  and  electrical  conductivity.  These  simpli¬ 
fications  enable  a  drastic  improvement  in  computational  efficiency 
while  sufficiently  capturing  the  on-cell  spatial  distribution  of  tem¬ 
perature  and  the  contribution  to  ohmic  loss  between  cells. 

With  these  assumptions,  we  derive  a  model  where  property  val¬ 
ues  along  a  ID  path  along  the  flow  direction  are  calculated  for 
four  domains  in  each  cell:  the  fuel  stream,  the  air  stream,  the 
anode-electrolyte-cathode  structure,  and  the  interconnect.  The 
quasi-2D  model  is  then  assembled  by  stacking  multiple  ID  paths 
through  the  thickness  direction  to  simulate  a  multi-cell  SOFC  stack. 
The  associated  computational  grid  for  the  cell  is  schematically 
depicted  in  Fig.  5(a).  This  arrangement  is  best  characterized  as  a 
quasi-2D  solution,  but  the  model  will  subsequently  be  referred  to 
here  as  SOFC-MP  2D  (to  differentiate  it  from  the  existing  3D  model). 

In  the  electrochemistry  model,  the  PEN  layers  are  represented 
with  one  computational  point  where  a  current-voltage  relation¬ 
ship  is  implemented.  Due  to  the  low  ionic  conductivity  of  the 
electrolyte  relative  to  the  high  electrical  conductivity  of  the  anode 
and  interconnect,  the  ohmic  voltage  drop  in  the  electrolyte  is  much 
greater  than  the  ohmic  drop  in  the  anode  and  interconnect.  This 
permits  the  code  to  assume  that  an  entire  planar  PEN  is  at  the  same 
working  voltage.  The  working  voltage  for  each  cell  in  the  stack  can 
then  be  calculated  without  the  need  for  a  full  electric  field  computa¬ 
tion.  The  current-voltage  relationship  at  each  computational  point 
is  then  used  to  obtain  the  cell’s  current  density  distribution  based 
on  the  local  species  concentrations,  temperature,  and  cell  voltage. 

Finally,  we  assume  that  the  stack  operates  at  steady  state  with 
the  initial  configuration  and  virgin  material  properties.  In  other 
words,  the  model  will  solve  the  electrochemistry  problem  with¬ 
out  considering  any  structural  changes  or  material  degradations 
caused  by  the  high  operating  temperature  of  the  stack.  The  ther¬ 
mal  deformations/stresses  at  the  component  level  and  degradation 
mechanism  processes  at  the  microstructural  level  in  the  fuel  cell  are 
other  important  areas  pursued  by  SOFC  researchers  because  the 
structural  response  obviously  depends  on  the  electrochemical  per¬ 
formance.  Usually  though,  structural  deformations  are  assumed  to 
be  small  enough  such  that  the  effect  on  gas  flow  and  electrochem¬ 
ical  performance  is  minimal. 

3.  Governing  equations 


stack/cell  level  models  (e.g.,  [13,44]).  Overall,  the  electrochem¬ 
istry  (EC)  model  provides  the  user  with  the  power  output,  species 
utilizations,  and  temperature  results  that  are  necessary  to  make 
engineering  decisions  on  the  design  and  operation  of  multi-cell 
stacks. 

3.1.  Electrochemistry 

The  EC  model  in  SOFC-MP  simulates  the  cell  electrochemical 
behavior  through  a  current-voltage  relationship,  which  provides 
the  cell  operating  voltage  as  a  function  of  current  density,  temper¬ 
ature,  cell  properties,  and  local  species  concentrations.  The  general 
form  of  the  current-voltage  relationship  consists  of  the  Nernst  volt- 
age  (Wiemst)  minus  the  ohmic  ('/ohmic),  activation  (^activation),  and 
concentration  (^concentration)  polarizations: 


Vcell  —  V^femst  ?7ohmic  ^activation  ^concentration  — /(Q 


(1) 


The  Nernst  voltage  for  a  hydrogen  fuel  consists  of  the  reversible 
open  circuit  voltage,  computed  from  the  change  of  Gibbs  free 
energy  for  H2  oxidation  at  standard  pressure  at  the  operating  tem¬ 
perature,  with  adjustment  for  the  actual  system  pressure  and  gas 
concentrations  used.  The  Nernst  potential  assuming  atmospheric 
fuel  and  oxidant  pressures  is  given  by 


^  +  ^1" 


/Ph2Po21/2\ 
V  Ph2  o  ) 


(2) 


where  A  Gf  is  the  enthalpy  change  of  formation  for  the  reaction  at 
the  local  temperature,  F  is  the  Faraday  constant,  RG  is  the  gas  con¬ 
stant,  T  is  the  temperature,  and  pH2,  Po2,  Ph2o  are  the  hydrogen, 
oxygen,  and  water  partial  pressures,  respectively.  Two  methods 
for  inputting  the  current-voltage  ( I-V )  relationship  are  imple¬ 
mented.  First,  the  approach  used  previously  in  other  modeling 
activities  [29,45]  is  implemented  as  the  default  calculation  method. 
This  form  shown  in  Eq.  (3)  calculates  the  cell  operating  voltage 
by  computing  terms  corresponding  to  the  local  Nernst  potential, 
ohmic  resistance,  activation  polarization,  cathode  02  concentra¬ 
tion  polarization,  anode  H2  concentration  polarization,  and  anode 
IT20  concentration  polarization  as 


The  governing  equations  describing  the  electrochemistry,  heat 
transfer,  and  fluid  flow  are  given  below  along  with  a  brief  descrip¬ 
tion  of  the  underlying  assumptions.  More  in-depth  discussion  on 
these  topics  can  be  found  in  the  references  on  system  models  and 


where  VNemst  is  the  Nernst  potential,  i  is  the  current  density,  is  the 
Ohmic  resistance,  bsinh_1(i/2l0)  is  the  activation  polarization,  RG 
is  the  gas  constant,  T  is  the  temperature,  F  is  the  Faraday  constant, 
io2  is  the  oxygen  transport  limiting  current,  i'h2  is  the  hydrogen 
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transport  limiting  current,  and  p^Q  and  p°2  are  the  water  and 
hydrogen  partial  pressures,  respectively.  The  activation  polariza¬ 
tion  is  described  by  a  Butler-Volmer  formulation  where 


b 


RGT 

aF 


(4) 


where  a  is  the  transfer  coefficient,  and  the  exchange  current  density 
i'o  is  given  by 


with  terms  that  provide,  respectively,  the  change  of  concentra¬ 
tion  and  density  p  with  time  t,  advection  due  to  the  velocity 
field  V ,  diffusive  flux  due  to  concentration  gradients},  and  species 
sinks/sources  Sz  The  EC  module  simplifies  the  gas  transport  prob¬ 
lem  by  assuming  the  gas  channels  are  incompressible,  steady-state, 
ID,  constant  velocity  flows  where  mass  diffusion  is  small  relative 
to  the  bulk  advection,  resulting  in  the  following  relation  for  species 
conservation: 


'0  =  P*eXP  (r|) 


(5)  PV  -k3. 


Ac,- 

pvA^rSi 


(8) 


where  Ea  is  the  activation  energy  and  Px  is  a  pre-exponential 
constant.  The  limiting  currents  calculated  here  assume  that  the 
dominant  mode  of  gas  transport  in  the  anode  and  the  cathode  is 
by  bulk  diffusion  through  uniform  microstructures  of  character¬ 
istic  porosity  and  tortuosity  [13,44].  The  limiting  current  for  the 
specie  ix  is  given  by 


;  2  FDfPx 
x  RTtx 


(6) 


where  DJff  is  the  effective  diffusion  coefficient,  px  is  the  partial 
pressure,  tx  is  the  electrode  thickness,  and  subscript  x  denotes  the 
corresponding  species  of  H2,  02,  or  H20.  Finally,  the  resistance 
of  the  stack  structure  is  computed  from  the  layer  resistances  as 
Rt  =  Rp  +  Re  +  Rn+Rb ,  where  Rp  is  the  cathode  electrical  resistance,  Re 
is  the  electrolyte  ionic  resistance,  Rn  is  the  anode  electrical  resis¬ 
tance,  and  Rb  is  the  bipolar  plate  electrical  resistance. 

Since  Eq.  (3)  may  not  be  suitable  for  all  users  who  want 
to  implement  their  own  proprietary  electrochemical  routines  or 
experimental  data,  the  second  option  is  a  user-defined  electro¬ 
chemistry  subroutine.  The  subroutine  is  written  in  an  interpreted 
language  Lua  [49],  which  permits  the  subroutine  to  be  compiled 
directly  upon  execution  of  SOFC-MP.  This  permits  the  user- 
subroutine  to  be  of  any  mathematical  functional  or  tabular  form 
that  returns  only  the  appropriate  voltage  based  on  the  local  tem¬ 
perature,  species  partial  pressures,  and  current  density  arguments 
sent  from  the  main  solver  algorithm.  Details  of  the  Lua  interface  for 
user-defined  inputs  are  discussed  in  Section  3.3.2. 

The  EC  model  in  SOFC-MP  is  also  capable  of  simulating  elec¬ 
trochemical  reactions  with  composite  fuel  mixtures  containing  H2, 
H20,  CO,  C02,  N2,  and/or  CH4.  For  mixed  fuels,  it  is  assumed  in 
this  model  that  the  amount  of  CO  or  CH4  that  is  electrochemi- 
cally  oxidized  is  small  and  only  H2  is  oxidized  [44].  Li  and  Chyu 
[46]  showed  that  electrochemical  oxidation  of  both  CO  and  H2  at 
the  anode  yields  the  same  Nernst  potential  in  the  SOFC  as  long  as 
the  water-gas  shift  reaction  is  in  equilibrium.  This  was  confirmed 
experimentally  by  Weber  et  al.  [47]  who  observed  similar  electro¬ 
chemical  performance  with  direct  H2,  CO,  and  CH4  oxidation  as 
long  as  carbon  deposition  was  suppressed.  Matsuzaki  and  Yasuda 
[48]  showed  that  experimental  oxidation  rates  for  H2  were  two  to 
three  times  faster  than  CO  for  the  Ni:YSZ  anode/YSZ  electrolyte  sys¬ 
tem.  Therefore,  the  model  assumes  that  the  CH4  reformation  rate  is 
kinetically  controlled,  the  water-gas  shift  reaction  is  fast  and  always 
in  chemical  equilibrium,  and  that  fuel  cell  only  oxidizes  H2.  Since 
the  mixed  fuels  of  interest  in  this  work  are  typically  compositions 
with  a  significant  molar  fraction  of  H2  and  a  steam-to-carbon  ratio 
of  2-3  is  used  for  any  carbon  species  fuels,  the  same  formulation 
for  the  I-V  relationship  in  Eq.  (2)  can  thereby  be  used. 


3.2.  Conservation  of  mass/species 

The  generic  transport  equation  for  conservation  of  species  con¬ 
centration  Cj  is  given  by 

-gp  =  V  •  (pVq)  =  V  •  J  +  Sj  (7) 


where  p  is  the  density,  v  is  the  velocity  along  the  flow  direction,  and 
Si  is  the  species  sink/source  due  to  the  electrochemical  reactions. 
This  formulation  circumvents  the  need  to  solve  the  complete  fluid 
flow  problem  involving  the  Navier-Stokes  equations  and  hence 
results  in  significantly  improved  computational  efficiency. 

When  an  appropriate  I-V  model  is  specified,  the  current  density, 
i,  can  be  calculated  from  the  assumed  cell  voltage,  fuel  tempera¬ 
ture,  and  gas  species  composition.  The  overall  cell  reaction  for  the 
oxidation  of  hydrogen  is 

H2  +  io2  -*  H20  (9) 

The  corresponding  oxygen  and  hydrogen  consumption  rates  [50] 
for  a  stack  with  N  cells  are 

I  iNA  _i  .  I  iNA  _ i 

n°2  =  4 p  =  ~4F  moles  s  nu2  =  2 F=  2F  moles  s  (10) 

where  /  is  the  current  generated  for  the  specific  control  volume 
and  A  is  the  active  area.  Each  mole  of  oxygen  oxidizes  two  moles 
of  hydrogen  to  produce  two  moles  of  water  and  two  electrons.  The 
change  to  the  local  molar  flow  rates  for  the  hydrogen,  water,  and 
oxygen  species  can  then  be  easily  obtained. 

When  CO  is  added  to  the  fuel,  the  water-gas  shift  reaction  con¬ 
verts  CO  into  H2  according  to  the  following  reaction: 

CO  +  H20  H2  +  C02  (11) 


The  reaction  is  assumed  to  be  fast,  such  that  the  species 
composition  quickly  reaches  local  equilibrium  based  on  a 
temperature-dependent  equilibrium  constant  computed  from  the 
free  energy  of  the  reaction: 


Kshift,CO  = 


Ph2^C02 

Ph2oPco 


=  exp 


^  AGghift  ^ 


=  exp  - 


Eproducts^-^O-En 


M  -  TSi ) 


RT 


(12) 


where  the  change  in  free  energy  AGshift  is  computed  using  the 
changes  in  enthalpy  Hi  and  entropy  Sz  for  each  of  the  participating 
species  in  the  reaction.  With  both  hydrogen  and  carbon  monoxide 
in  the  fuel,  the  fuel  composition  along  the  flow  direction  is  obtained 
by  simultaneously  solving  the  shift  equilibrium  and  oxygen  con¬ 
sumption  rate  equation  for  each  control  volume  marching  along 
the  flow  direction.  Since  the  user  may  implement  any  arbitrary  fuel 
composition  into  the  model,  highly  non-equilibrated  fuels  contain¬ 
ing  CO  would  shift  rapidly  at  the  fuel  inlet  region  with  high  heat 
outputs.  The  model  also  includes  a  user-defined  distance  param¬ 
eter  over  which  the  fuel  composition  is  allowed  to  progress  from 
the  initial  state  to  full  equilibrium  where  the  local  value  of  Kshift  co 
changes  linearly  from  the  value  at  the  inlet  to  the  equilibrium  value. 
Generally,  testing  with  various  realistic  fuel  compositions  showed 
few  numerical  issues  due  to  other  solution  control  mechanisms 
that  aid  stable  convergence.  If  methane  is  added  to  the  fuel,  steam 
reformation  occurs,  which  is  an  endothermic  reaction  where  CH4 
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is  catalyzed  by  the  nickel  in  the  anode  cermet: 

CH4  +  H20^3H2  +  C0  (13) 

The  reformation  reaction  is  kinetically  controlled.  Rate  equa¬ 
tions  are  difficult  to  obtain  experimentally,  but  various  relation¬ 
ships  for  Ni:YSZ  anodes  have  been  established  from  experiments 
for  use  in  modeling  efforts.  Xu  and  Froment  [51],  Achenbach  and 
Riensche  [52,42],  Dicks  et  al.  [53],  and  Ahmed  and  Foger  [54]  cal¬ 
culated  the  reaction  rate  of  methane  according  to  an  empirical 
formula  derived  from  experiments.  At  cell  operating  temperatures, 
the  methane  reforming  occurs  relatively  quickly  until  the  composi¬ 
tion  reaches  equilibrium.  The  Achenbach  and  Riensche  relation  [  52  ] 
is  implemented  as  the  default  reforming  rate  equation,  but  any  user 
relation  can  be  implemented  through  a  Lua  subroutine.  The  equi¬ 
librium  constant  is  defined  similarly  to  that  for  the  water-gas  shift 
reaction  as 


^shift,CH4  = 


Ph2Pco 

Pch4Ph2o 


(14) 


and  the  value  of  /CShift,cH4  is  calculated  from  the  change  in  free 
energy  for  the  reaction  components  in  the  same  way  as  described 
above  for  CO. 

SOFC-MP  provides  three  approaches  for  the  simulation  of 
methane  reforming.  Because  methane  reforming  is  an  endother¬ 
mic  reaction,  the  temperature  field  and  simulation  result  near  the 
fuel  inlet  area  are  sensitive  to  the  reforming  model  chosen.  The 
first  approach  applies  a  default  reforming  rate  as  described  above, 
e.g.,  one  provided  by  Achenbach  [42].  The  methane  reforms  along 
the  flow  path  according  to  the  computed  reaction  rate  until  equilib¬ 
rium  is  reached,  and  the  gas  compositions  remain  fully  equilibrated 
for  the  rest  of  the  fuel  flow.  In  addition  to  the  default  reforming 
rate,  the  model  also  allows  users  to  define  their  own  reforming 
rate  expression  through  a  Lua  interface,  which  will  be  discussed 
in  Section  3.3.2.  Another  approach  is  to  assume  that  the  methane 
reforming  reaction  is  also  relatively  fast  and  can  be  approximated 
with  the  equilibrium-path-length  approximation  described  above 
for  CO.  For  example,  the  fuel  composition  described  in  Section  5.1 .2 
that  initially  contains  1 7%  methane  has  an  equilibrium  value  of  less 
than  0.01%  methane  at  high  operating  temperatures  above  900  °C. 
The  numerical  result  from  the  reforming  rate  approach  using  the 
Achenbach  rate  expression  indicates  that  the  composition  of  CFI4 
is  90%  consumed  at  about  5-10%  of  the  10  cm  cell  length,  although 
full  equilibrium  is  not  reached  until  20%  of  cell  length.  These  results 
indicate  that  the  equilibrium  path  approach  with  a  reasonable  per¬ 
centage  cell  length  based  on  available  experimental  observations 
could  be  a  good  alternative  when  a  rate  expression  is  difficult  to 
obtain.  By  balancing  the  fuel  species  to  equilibrium  according  to 
the  current  density  calculated  by  the  l-V  relationship,  the  coupling 
between  the  fluid  flow  and  the  electrochemical  reactions  is  thus 
achieved,  and  it  is  used  to  satisfy  the  mass  and  species  balance  for 
each  computational  point. 


provided  by  the  fuel  and  oxidant  flows  and/or  exterior  losses  to 
the  environment.  Typically  in  most  designs,  oxidant  flow  rates  are 
increased  to  provide  sufficient  cooling  through  forced  convection. 
The  heat  transfer  mechanisms  included  in  the  model  are  (1)  con¬ 
duction  in  the  cell  and  interconnect,  (2)  convection  between  the 
cell/interconnect  and  the  fuel/oxidant  flows,  (3)  advection  within 
the  fuel/oxidant  flows,  (4)  radiation  between  the  cell  and  inter¬ 
connect,  (5)  radiation/convection  between  the  cell  and  insulated 
enclosure,  and  6)  radiation/convection  from  the  insulated  enclo¬ 
sure  to  the  external  environment. 

SOFC-MP  contains  a  core  algorithm  that  solves  the  steady-state 
heat  transfer  problem  by  the  finite-volume  method  to  provide 
the  temperature  distribution  in  the  entire  fuel  cell  stack.  Consid¬ 
ering  the  entire  model  domain,  the  enthalpy  changes  within  the 
stack,  the  electrical  power  output,  and  the  heat  loss  to  the  envi¬ 
ronment  are  balanced  until  the  steady-state  solution  is  achieved. 
The  EC  model  is  coupled  to  the  heat  transfer  model  through  the 
assigned  current-voltage  relationship  that  defines  the  usable  elec¬ 
trical  power  based  on  the  voltage  of  each  volume  representing 
the  cell.  The  enthalpy  of  the  gases  is  computed  using  the  standard 
Chemical  Equilibrium  with  Applications  (CEAs)  FORTRAN  property 
library  [56,57].  For  the  assigned  l-V  performance,  all  the  volumetric 
heat  generated  (or  absorbed)  from  the  H2  oxidation,  CO  shift,  and 
CFI4  reforming  reactions  as  well  as  the  Joule  heating  and  polariza¬ 
tion  losses  are  accounted  for  in  the  enthalpy  change  calculation. 

The  heating  rate  due  to  the  enthalpy  change  with  time,  ( AH/ A  t), 
from  advection  of  species  through  a  computational  volume  or  the 
ionic  transport  of  the  02  is  then 

^  =  (is) 


where  hj  is  a  molar  flow  rate,  Mj  is  the  molar  mass,  hj  is  the  specific 
enthalpy,  and  j  is  summed  over  the  number  of  present  species.  The 
heat  transfer  rates  between  adjacent  volume  elements  are  com¬ 
puted  with  conventional  equations  for  solid  conduction  Eq.  (16), 
solid-gas  convection  Eq.  (17),  solid-solid  gap  radiation  Eq.  (18),  and 
radiation  to  the  ambient  Eq.  (19): 


^conduction  — 


kA 


T-Tn 

Ax 


(16) 


*7 convection  —  hcA(T  Tn )  (17) 

*7gap_radiation  —  ~AS\  £2<tAT^(T  —  Tn)  (18) 


9ambient_radiation  —  4scrAT  ( T  Ta ) 


(19) 


where  k  is  thermal  conductivity,  hc  is  the  convective  film  coefficient, 
T  is  the  local  temperature,  Tn  is  the  temperature  of  the  neighbor¬ 
ing  volume  element,  Ta  is  the  ambient  environment  temperature, 
A  is  the  area  of  the  element  face  through  which  the  heat  flux  is 
calculated,  Ax  is  the  conduction  length,  e  is  emissivity,  and  a  is 
the  Stefan-Boltzmann  constant.  The  usable  electrical  power  Pelec  is 
simply 


3.3.  Conservation  of  energy 


^elec  —  17iAcen 


(20) 


The  temperature  distribution  of  an  SOFC  has  a  very  important 
effect  on  its  performance  (e.g.,  [50,55])  since  the  thermal-physical 
properties,  electrochemistry  performance,  and  gas  species  equi¬ 
librium  composition  is  all  tightly  coupled  with  the  temperature. 
Solving  the  thermal  energy  balance  is  thus  one  of  the  major  tasks 
in  the  SOFC-MP  model.  Fleat  is  primarily  generated  from  the  elec¬ 
trochemistry  reactions  where  oxidation  of  H2  and  the  CO  water-gas 
shift  are  exothermic  reactions  while  CH4  steam  reformation  is 
endothermic.  Joule  heating  is  also  created  in  the  solids  because  of 
ohmic  resistance.  To  prevent  excessive  temperatures  in  the  stack, 
this  generated  heat  must  be  dissipated  through  convective  cooling 


where  V  is  the  cell  operating  voltage,  i  is  the  cell  average  current 
density,  and  Acell  is  the  cell  active  area.  Therefore,  the  energy  bal¬ 
ance  for  the  generic  finite  volume  is  thereby  established  by 


+  ^elec  —  6 


(21) 


where  only  the  terms  corresponding  to  the  appropriate  thermal 
transport  mechanisms  (Eq.  (15)  through  Eq.  (20))  active  in  each 
volume  of  interest  are  used,  and  energy  is  locally  conserved  in 
every  volume  element.  The  electrochemical  reactions  are  assumed 
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to  occur  within  the  anode,  so  the  heat  from  the  enthalpy  change  is 
assigned  to  the  PEN  structure. 

During  the  solution  process,  the  maximum  temperature  differ¬ 
ence  between  iterations  for  every  point  is  checked  until  it  is  reduced 
below  a  user-defined  convergence  tolerance.  The  global  energy  bal¬ 
ance  for  the  entire  stack  can  then  be  validated  after  the  solution 
by  comparing  the  net  chemical  energy  Qreaction.  useful  electrical 
power  Pelec,  enthalpy  changes  of  the  outlet  fuel/oxidant  streams 
Qgases,  and  exterior  heat  losses  Qi0Ss  according  to  the  following 
relation: 

AH 

A^imbalance  =  ^  +  Pelec  +  Qloss  —  Qreaction  +  Qgases 

+  Peiec  +  Qloss  ^  0  (22) 

3.3 A.  Boundary  conditions 

The  model  computes  the  stack  temperature  distribution  based 
on  the  input  properties  of  the  fuel  and  air  streams,  the  operating 
voltage,  and  the  heat  loss  to  its  surroundings.  The  inlet  composition, 
inlet  temperature,  and  flow  rates  of  the  fuel  and  air  streams  are 
defined  by  the  user. 

The  model  has  the  capability  to  define  external  thermal  bound¬ 
ary  conditions  that  mimic  the  cell  placement  within  a  furnace 
or  an  insulated  enclosure.  These  boundary  conditions  account 
for  heat  transfer  by  convection  and  radiation  from  the  stack  to 
its  surroundings.  For  the  furnace  condition,  the  convective  heat 
transfer  coefficient,  ambient  air  temperature,  stack  emissivity,  and 
furnace  temperature  are  defined  for  the  exterior  of  the  stack.  Val¬ 
ues  may  be  assigned  independently  for  the  bottom,  top,  front, 
back,  and  sides  of  the  stack.  If  the  stack  is  contained  within 
an  insulated  enclosure,  additional  values  for  the  insulation  ther¬ 
mal  conductivity,  thickness,  convection  coefficient  between  the 
stack  and  insulation,  and  emissivity  of  the  insulation  are  also 
defined  for  each  side  in  addition  to  the  previous  external  bound¬ 
ary  conditions.  Furthermore,  variation  of  the  inlet  temperature 
of  air  and  fuel  from  the  bottom  to  the  top  of  the  stack  can 
be  specified  to  best  mimic  the  actual  operating  environment  if 
desired. 

3.3.2.  User-defined  inputs 

Different  cell  designers  will  likely  have  different  modeling 
approaches  for  simulating  the  I-V  curve  for  their  cells.  For  flex¬ 
ibility,  a  well  designed  ability  to  work  with  a  broad  spectrum 
of  user-defined  electrochemistry  models  is  desirable.  Therefore,  a 
general  SOFC  modeling  program  must  not  be  constrained  by  using 
only  one  specific  model  and  equation  set.  The  capability  to  include 
customized  electrochemistry  models  is  implemented  in  the  SOFC- 
MP  model.  This  flexibility  also  works  well  on  the  principle  that 
SOFC-MP  users  do  not  desire  to  share  their  proprietary  informa¬ 
tion  with  software  developers  and  other  users.  For  this  purpose, 
programming  language  Lua  [49]  was  chosen  as  the  designated  lan¬ 
guage  for  user  electrochemistry  subroutines  that  interface  with 
SOFC-MP.  The  choice  of  Lua  was  made  because  of  its  relatively  sim¬ 
ple  C  application  program  interface  (API),  popularity,  ease  of  use, 
and  availability  of  support.  For  SOFC-MP,  the  primary  Lua  interface 
is  for  a  user  model  to  calculate  voltage  based  on  the  current  density, 
temperature,  and  species  composition  parameters,  but  other  inter¬ 
faces  have  been  developed  as  needed  to  implement  customized 
thermal-physical  properties  or  control  the  methane  reformation 
rate. 

4.  Numerical  algorithms  and  implementation 

The  2D  EC  module  for  the  fuel  cell  was  programmed  in  Visual 
Studio  C++  2008  for  Windows  and  GNU  C++  for  Linux,  interfacing 
with  many  subroutines  written  in  FORTRAN  77.  The  program  is 


currently  in  standalone  format,  but  it  can  be  readily  converted  into 
a  library  or  a  dynamic-link  library  (DLL)  to  be  integrated  with  a 
more  comprehensive  SOFC  system-level  model. 

The  multi-cell  stack  SOFC  is  modeled  in  2D  control  volumes  as 
shown  in  Fig.  5.  Assuming  the  number  of  cells  is  Ncens,  and  the 
numerical  increment  in  the  flow  direction  is  Nx,  there  will  be  a 
total  of  (4Nceils  + 1  )NX  control  volumes  for  the  whole  SOFC  stack  in 
our  finite  volume  method  (FVM)  solution  scheme. 

The  solution  starts  with  an  initial  guess  on  the  temperature  field 
and  cell  voltage  based  on  inlet  air  and  fuel  temperature  as  well  as 
the  average  cell  voltage  specified  from  the  input  file.  It  then  involves 
four  numerical  steps  in  two  iteration  loops. 

The  model  follows  a  control  volume  marching  in  the  fuel  flow 
direction  from  the  control  volumes  in  the  inlet  area  toward  the 
outlet  area.  In  each  control  volume,  the  air  and  fuel  gas  composition 
is  known  from  the  previous  control  volume.  With  the  oxygen  and 
hydrogen  consumption  rates  calculated  according  to  Eq.  (10),  the 
remaining  species  content  will  be  calculated,  and  the  gas  partial 
pressure  will  be  balanced  to  give  the  new  fuel  and  air  composition 
for  the  subsequent  volume.  For  the  same  control  volume,  the  model 
uses  the  default  internal  relation  Eq.  (3)  or  another  user-specified 
I-V  relationship  to  calculate  the  current  density  based  on  the  gas 
and  oxidant  temperature,  species  composition,  and  cell  voltage. 

A  thermal  equilibrium  is  required  for  each  control  volume  as 
shown  in  Fig.  5(a).  Each  control  volume  typically  has  four  neigh¬ 
boring  control  volumes,  and  with  each  one,  heat  conduction, 
convection,  or  radiation  takes  place.  All  of  these  heat  transfer  modes 
involve  the  temperature  values  of  the  current  control  volume  and 
its  neighbors.  When  all  mechanisms  are  considered,  including  the 
heat  generated  from  the  electrochemistry  reaction  and  the  heat 
transfer  to  the  environment,  the  equation  to  be  solved  is  of  the 
form: 

[Ai\.[Ti\=Fi  (23) 

Here  [A/]  is  a  1  x  5  matrix  with  coefficients  for  the  control  volume 
and  its  neighbors  from  thermal  calculation,  [Tz  ]  is  a  5x1  matrix 
representing  the  temperature  at  the  center  of  each  control  volume, 
and  Fj  is  the  resulting  force  term  for  this  thermal  equilibrium  equa¬ 
tion.  When  equations  for  all  control  volumes  are  combined,  a  matrix 
equation  of  temperatures  is  formed: 

[A].[T]  =  [F]  (24) 

Here  [A]  is  a  square  matrix  with  the  size  of  the  number  of  control 
volumes,  [T]  is  the  temperature  vector  for  all  control  volumes,  and 
[F]  is  the  vector  of  the  force  term.  Because  both  [A]  and  [F]  depend 
on  the  temperature  field  for  all  control  volumes,  it  takes  multi¬ 
ple  iterations  to  achieve  a  solution.  The  convergence  is  considered 
achieved  when  the  temperature  difference  at  every  control  volume 
from  two  consecutive  iteration  steps  is  within  the  user  preset  tol¬ 
erance,  e.g.,  0.1  °C.  The  calculation  of  [T]  for  a  given  cell  voltage 
distribution  constitutes  the  inner  iteration  for  the  model. 

For  this  given  cell  voltage  distribution,  the  current  for  each  cell 
in  the  stack  is  calculated  by  integrating  the  current  density  along 
the  flow  direction.  Because  all  cells  are  connected  in  series  in  the 
same  circuit,  each  cell  must  have  the  same  total  current.  The  differ¬ 
ence  of  current  in  each  cell  indicates  that  the  cell  voltage  must  be 
redistributed.  This  voltage  rebalance  activity  forms  the  outer  loop 
iteration  that  converges  when  the  current  difference  among  cells  is 
within  the  user-specified  tolerance,  e.g.,  1%  of  the  average  current. 
The  model  reaches  a  sound  solution  when  convergences  for  both 
the  inner  temperature  loop  and  the  outer  current-voltage  loop  are 
achieved.  For  a  typical  tolerance  of  0.1  °C  for  temperature  and  1% 
for  current,  it  takes  fewer  than  20  iterations  on  the  inner  tempera¬ 
ture  loop  and  less  than  10  iterations  on  the  outer  current-voltage 
loop  to  converge.  The  number  of  iterations  needed  for  the  inner 
loop  then  decreases  as  the  outer  loop  converges.  The  number  of 
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iterations  needed  for  a  counter-flow  case  is  approximately  double 
that  of  a  co-flow  case. 

The  multi-physics  algorithm  was  constructed  on  a  given  average 
cell  voltage  on  a  multi-cell  stack.  Several  other  options  were  derived 
by  iteratively  applying  the  same  algorithm: 

•  Solve  the  multi-stack  model  based  on  a  given  average  current 
density. 

•  Solve  the  multi-stack  model  based  on  a  given  fuel  utilization  rate. 

•  Generate  a  stack  wise  I-V  curve  on  fixed  fuel  and  air  flow  rate. 

•  Generate  a  stack  wise  I-V  curve  on  fixed  fuel  and  air  utilization 
rate. 

For  all  computational  options,  the  following  input  information 
must  be  specified  by  the  user: 

•  Cell  geometry  parameters  (number  of  cells,  cell  length,  width,  and 
height;  thickness  of  interconnect,  PEN). 

•  Fuel  and  air  thermodynamic  state  (temperature  and  pressure) 
and  species  compositions  at  the  inlet  area. 

•  Boundary  conditions  (temperatures  of  the  surrounding  top, 
bottom,  front,  back,  and  side  enclosures  including  thermal  con¬ 
ductivity,  emissivity,  film  coefficients,  etc.). 

•  Operating  conditions  of  either  the  desired  average  working  volt¬ 
age  of  all  cells  or  the  desired  average  current  density  of  all  cells 
is  defined  based  on  the  designated  solution  mode. 

The  model  solves  the  multi-physics  system  and  predicts  the 
following  information  for  the  SOFC  system: 

•  Overall  SOFC  stack  performance:  total  power  output,  fuel  and  air 
utilization,  heat  loss  at  each  side  of  the  boundary,  working  voltage 
for  each  cell  in  the  stack,  working  current  density. 

•  Detailed  profiles  of  various  thermophysical  parameters:  temper¬ 
ature,  current  density,  species  composition,  heat  generation  at 
each  control  volume  (fuel,  air,  PEN,  or  interconnect). 

In  addition  to  the  input  parameters  related  to  the  operating  con¬ 
ditions,  users  also  have  the  flexibility  to  control  the  parameters  for 
numerical  iteration,  such  as  maximum  values  for  the  number  of 
inner  loop  (temperature)  and  outer  loop  (current-voltage)  itera¬ 
tions,  tolerances  for  both  iterations,  and  other  tuning  parameters 
such  as  relaxation  factors  for  both  temperature  and  current- voltage 
iterations.  Furthermore,  the  model  offers  a  restart  capability  by 
saving  temperature  and  current  values  from  the  previous  simu¬ 
lation.  This  restart  feature  has  proven  to  be  a  good  feature  that 
will  shorten  the  computation  time  when  small  changes  are  made 
to  input  parameters  (e.g.,  during  successive  calls  during  a  system 
simulation). 

4.1  A.  Results  and  post-processing 

Once  the  convergence  criteria  are  satisfied,  the  final  results  are 
tabulated  and  output  for  easy  access  by  the  user.  Distributions  for 
the  current  density,  temperature,  and  species  concentrations  along 
the  path  length  and  voltage  through  the  stack  height  are  output  for 
each  cell  in  the  stack.  The  energy  balance  is  output  to  indicate  the 
amount  of  total  chemical  energy  that  is  converted  to  useful  power, 
expended  as  enthalpy  increases  in  the  fuel/oxidant  flows,  or  lost  to 
the  external  environment.  The  energy  balance  also  assures  the  user 
that  convergence  has  been  achieved  with  a  small  tolerable  energy 
imbalance.  Derived  stack  performance,  such  as  power  output,  fuel 
utilization,  air  utilization,  maximum  cell  temperature,  cell  temper¬ 
ature  difference,  gas  outlet  temperatures,  and  maximum  current 
density  are  output  for  the  user.  From  these  sets  of  data  and  distri¬ 
bution  results,  the  user  can  calculate  any  additional  performance 


metrics  of  interest  for  the  stack.  These  post-processing  computa¬ 
tions  have  been  implemented  into  an  Excel  dashboard  worksheet 
to  automatically  calculate  all  the  stack  engineering  metrics  and 
plot  all  the  distributions  of  results  using  macros.  This  provides  a 
very  simple  procedure  for  the  user  to  perform  stack-level  engi¬ 
neering  calculations  to  satisfy  the  design  constraints  and  associated 
operation  criteria. 


5.  Model  verification  and  validation 

A  standard  set  of  benchmark  cases  for  SOFC  operation  and 
modeling  have  not  been  formally  adopted  by  the  fuel  cell  com¬ 
munity,  but  some  benchmark  cases  previously  established  through 
collaborative  activities  at  the  International  Energy  Agency  (IEA) 
have  emerged  as  cases  of  interest.  Modeling  predictions  from 
a  nine-member  round-robin  test  were  compared  by  Achenbach 
[42]  for  co-flow,  counter-flow,  and  cross-flow  geometries  with 
non-reforming  and  reforming  fuels.  Subsequently,  these  bench¬ 
mark  cases  have  informally  served  as  test  cases  for  various  other 
researchers  who  have  reported  results  for  their  SOFC  simulation 
tools  [28,43,59-61  ].  In  contrast  to  the  current  state-of-the-art  high 
performance  anode-supported  SOFCs,  the  benchmark  simulates  an 
electrolyte-supported  cell  with  ceramic  interconnect  and  a  low  cur¬ 
rent  density  of  3000 Am-2.  Regardless,  the  results  can  still  serve 
as  valuable  validations  for  the  modeling  procedures  developed  in 
this  study.  Therefore,  since  high  fidelity  experimental  data  from 
well  defined  stacks  has  not  yet  been  obtained  for  public  use  due  to 
proprietary  restrictions,  the  model  predictions  from  this  study  are 
compared  to  other  modeling  results  in  the  literature  for  validation 
of  the  approach. 

The  modeling  efforts  described  above,  however,  does  not  seem 
to  have  companion  cases  for  multi-cell  stacks.  Rather,  it  appears 
that  few  multi-cell  stack  models  exist  in  the  literature  that  fully 
describe  the  model  basis,  inputs,  material  properties,  and  results 
output  with  sufficient  technical  detail  for  replication.  One  example 
which  is  adequately  described  and  will  be  used  here  is  the  five-cell 
stack  studies  reported  by  Burt  et  al.  [62].  Their  modeling  approach 
is  similar  to  that  taken  here  and  consists  of  a  series  of  stacked  ID 
solutions  to  obtain  the  2D  variation  through  the  stack  cross-section. 
Their  examples  evaluated  both  anode  and  electrolyte  supported 
cells  with  parametric  studies  on  the  effects  of  flow  maldistribution 
and  radiation  heat  transfer.  Flere,  four  cases  using  five-cell  stacks 
and  uniform  flow  distributions  will  be  simulated.  The  anode  and 
electrolyte  supported  cells  will  be  considered  with  and  without 
radiation  between  the  cell  and  interconnect,  and  comparisons  of 
the  predicted  temperature  profiles  will  be  made. 

5.1.  Single  cell  benchmark  model 

For  the  IEA  cases,  the  baseline  conditions  for  the  model  include 
3000  A  nrr2  current  density,  900  °C  inlet  gas  temperatures,  85%  fuel 
utilization,  a  stoichiometric  air  ratio  of  7.0,  and  adiabatic  boundary 
conditions.  The  IEA  electrochemistry  model  assumes  that  the  anode 
and  cathode  concentration  polarizations  are  equal  to  the  electrolyte 
ohmic  polarization  for  the  sake  of  simplicity.  Our  quasi-2D  model 
can  be  compared  with  co-  and  counter-flow  geometries  for  each 
case.  For  each  case,  the  fuel  gas  composition  is  entered  as  a  known 
parameter,  and  the  IEA  electrochemistry  model  is  replicated  and 
implemented  through  the  Lua  program  interface.  The  benchmark 
cases  are  adiabatic,  meaning  that  no  heat  is  lost  to  the  environ¬ 
ment,  but  the  convection  coefficient  and  radiation  emissivity  for 
the  solid-fluid  interfaces  are  estimated.  Output  metrics  for  compar¬ 
ison  include  voltage,  power,  efficiency,  current  density  distribution, 
temperature  distribution,  and  air/fuel  outlet  temperatures. 
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Fig.  6.  Simulation  results  for  (a-c)  Case  1  co-flow  benchmark  and  (d-f)  Case  2  counter-flow  benchmark  without  methane  reforming. 


5.1  .1.  Cases  1  and  2:  humidified  hydrogen  fuel 

For  Case  1  (co-flow)  and  Case  2  (counter-flow),  the  fuel  is  90% 
H2  while  the  oxidant  is  standard  air  with  21%  02,  both  with  inlet 
temperatures  of  900  °C.  The  total  current  density  is  assigned  in  the 
model,  and  the  flow  rates  are  adjusted  to  achieve  the  correct  uti¬ 
lizations.  The  resulting  distributions  for  the  current  density,  cell 
temperature,  and  fuel  composition  for  the  co-flow  and  counter¬ 
flow  geometries  are  shown  in  Fig.  6.  In  the  co-flow  geometry, 
the  current  density  is  relatively  uniform  over  the  first  half  of  the 
cell.  This  is  due  to  the  competing  factors  along  the  cell  length  of 
improved  electrochemical  activity  due  to  lower  ohmic  losses  as  the 
temperature  gets  higher  versus  the  reduced  Nernst  potential  with 
reduced  H2  concentration  due  to  consumption.  Near  the  outlet, 
the  temperature  is  high,  but  the  fuel  composition  is  low,  result¬ 
ing  in  lower  current  density.  For  the  counter-flow  case,  the  highest 
cell  temperature  and  the  highest  FI2  concentration  result  in  much 
higher  current  density  at  the  fuel  inlet  region.  The  temperature  is 
highest  here  because  of  the  hot  oxidant  stream,  which  removes 
most  of  the  electrochemical  heat  from  the  cell.  The  current  density 
then  drops  dramatically  along  the  cell  to  reach  a  minimum  at  the 


air  inlet  side  where  both  cell  temperature  and  H2  concentration  are 
lowest. 

5.1.2.  Cases  3  and  4:  reforming  fuel 

For  Case  3  (co-flow)  and  Case  4  (counter-flow),  the  fuel  composi¬ 
tion  is  based  on  a  50%  pre-reformed  methane  composition  while  the 
oxidant  is  again  standard  air.  The  resulting  input  molar  composition 
for  the  fuel  is  then  26.26%  H2, 49.34%  H20, 2.94%  CO,  4.36%  C02,  and 
17.10%  CFI4.  In  our  model,  the  reforming  rate  model  established  by 
Achenbach  and  Riensche  [52]  is  used  to  calculate  the  composition  of 
CH4  and  other  species  until  full  equilibrium  is  reached.  The  results 
for  the  co-flow  and  counter-flow  geometries  are  shown  in  Fig.  7.  For 
the  co-flow  case,  the  current  density  is  fairly  uniform  across  the  cell. 
Due  to  the  rapid  reformation  near  the  fuel  inlet  region,  the  FI2  con¬ 
centration  is  high,  but  the  temperature  decreases  significantly  from 
the  endothermic  reaction,  resulting  in  slightly  reduced  current  pro¬ 
duction.  The  current  density  increases  as  the  cell  becomes  hotter 
downstream.  Near  the  fuel  outlet,  the  low  fuel  concentration  results 
in  a  low  Nernst  value  even  though  the  temperature  is  highest.  For 
counter-flow,  the  methane  still  gets  reformed  quickly  near  the  fuel 
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Fig.  7.  Simulation  results  for  (a-c)  Case  3  co-flow  benchmark  and  (d-f)  Case  4  counter-flow  benchmark  with  methane  reforming. 


inlet  side,  but  the  heated  air  stream  compensates  for  the  reforma¬ 
tion  endotherm,  resulting  in  a  higher  current  density  peak  from  the 
high  temperature  and  high  H2  concentration.  This  is  similar  to  the 
case  without  reformation,  except  that  the  endotherm  pushes  the 
location  of  maximum  current  density  downstream  slightly. 

5. 1.3.  Single  cell  benchmark  comparison 

The  results  from  the  four  cases  using  the  present  model  are 
shown  in  comparison  to  other  results  from  the  literature  in  Fig.  8 
[42,60].  The  extreme  values  for  these  results  are  compared  because 
the  actual  distributions  for  most  of  these  metrics  are  unavailable, 
but  the  comparisons  made  to  available  distributions  such  as  Achen- 
bach’s  temperature  data  [43  ]  shown  in  Fig.  9  have  been  good  and  the 
method  of  comparison  is  deemed  to  be  acceptable.  Overall,  these 
results  are  generally  within  the  ranges  of  the  other  models  and  are 
thereby  considered  to  be  an  acceptable  match,  especially  with  the 
approximations  utilized  to  simulate  the  rib  geometry.  The  electri¬ 
cal  performance  results  were  very  comparable.  The  model  voltages 
were  close  to  the  mean  value  of  the  reference  predictions,  while  the 
minimum  and  maximum  current  densities  fell  within  the  range  of 
the  literature  values.  This  suggests  that  the  overall  current  density 


distributions  and  fuel  distributions  across  the  cell  are  suitably  cap¬ 
tured  by  this  model.  The  temperature  results  compared  reasonably 
well  to  the  reference  results  also,  although  there  were  some  small 
discrepancies.  The  maximum  cell  temperatures  with  the  reform¬ 
ing  fuel  were  slightly  lower  than  the  reference  range  by  4-17  °C 
for  the  co-flow  configuration  (Case  3)  and  3-30  °C  for  the  counter¬ 
flow  configuration  (Case  4).  Also,  the  fuel  outlet  temperature  for 
the  co-flow  configuration  (Case  3)  was  4-9  °C  lower  than  the  ref¬ 
erence  range  and  the  air  outlet  temperature  for  the  counter-flow 
configuration  (Case  4)  was  2-12  °C  lower  than  the  reference  range. 
The  exact  reasons  for  these  discrepancies  remain  difficult  to  diag¬ 
nose,  but  can  possibly  be  attributed  to  the  following  factors:  (1 )  the 
shape  of  the  current  density  profile  determined  with  the  reform¬ 
ing  rate  equation,  (2)  the  convective  transport  properties,  and/or 
(3)  the  model  mesh  density.  The  rapid  methane  reformation  at  the 
fuel  inlet  side  results  in  steep  gradients  of  the  hydrogen  species 
concentration.  The  coupled  current  density  and  thermal  profile 
will  then  depend  on  the  heat  removed  by  the  thermal  transport 
mechanisms.  Since  the  assumptions  and  temperature  dependence 
for  the  convective  film  coefficients  are  also  unknown,  the  convec¬ 
tive  heat  transfer  to  the  gas  streams  may  be  slightly  different  from 
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Fig.  8.  Simulation  results  for  (a)  minimum/maximum  current  density,  (b)  minimum/maximum  cell  temperature,  (c)  voltage,  (d)  fuel  outlet  temperature,  and  (e)  air  outlet 
temperature  for  benchmark  Cases  1-4  with  comparisons  to  literature  results. 


the  reference  models.  Furthermore,  a  model  mesh  density  of  1000 
increments  was  used  for  these  benchmark  cases,  while  the  mesh 
resolution  of  the  reference  cases  is  unknown.  Mesh  resolution  dif¬ 
ferences  could  be  important  for  accurate  capture  of  peak  values 
near  the  edge  of  the  model  domain  such  as  fuel  and  air  outlet  tem¬ 
peratures.  In  Case  2  for  example,  the  fuel  temperature  changes 
significantly  from  the  900  °C  inlet  value  within  the  first  volume 
increment  (+16°C)  and  accordingly  affects  the  reformation  rate, 
while  in  Case  3  the  peak  cell  temperature  occurs  at  the  air  outlet 
edge.  Therefore,  despite  the  lack  of  the  necessary  detailed  assump¬ 
tions/inputs  to  replicate  the  original  benchmark  cases,  the  results 


still  compare  reasonably  well  overall  and  validate  the  use  of  this 
code  as  an  engineering  stack-level  simulation  tool. 

5.2.  Multi-cell  stack  example  model 

For  the  models  presented  in  Burt  et  al.  [62],  a  five-cell  stack 
was  used  consisting  of  900  cm2  co-flow  cell  with  the  10  cm  short 
dimension  along  the  flow  direction.  The  inlet  fuel  was  900  °C  wet 
hydrogen  with  3mol%  water,  while  the  inlet  oxidant  was  900  °C 
standard  air.  Current  levels  of  50-650  A  were  evaluated,  but  con¬ 
stant  fuel  and  air  flow  rates  were  used  so  that  utilizations  varied 
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Fig.  9.  Comparison  of  the  temperature  distribution  results  for  the  (a)  co-flow  and  (b)  counter-flow  cases  with  reforming  fuel. 


for  different  current  densities.  The  cell  was  assumed  to  be  adia¬ 
batic  with  no  heat  loss  to  the  external  environment.  The  relations 
for  the  electrochemistry  model  consisting  of  the  Nernst  voltage 
minus  ohmic,  activation,  and  concentration  polarizations  were 
coded  in  the  user-defined  Lua  interface.  Output  metrics  for  compar¬ 
ison  include  voltage,  power,  efficiency,  current  density  distribution, 
temperature  distribution,  and  air/fuel  outlet  temperatures.  In  this 
paper,  the  simulated  I-V  relationship  for  the  electrolyte-supported 
stack  and  the  temperature  distributions  for  the  top  and  bottom  cells 
of  the  five-cell  stack  are  used  for  comparisons. 

While  nearly  all  of  the  necessary  model  parameters  (e.g.,  dimen¬ 
sions,  thicknesses,  gap  heights,  flow  parameters,  EC  parameters, 
etc.)  for  comparison  were  available  in  the  associated  paper  and  Refs. 
[63-65],  detailed  analyses  found  three  parameters  that  required 
additional  assumptions  for  input  to  the  models.  First,  the  empir¬ 
ical  Nusselt  number  correlation  used  to  compute  the  convective 
coefficients  for  heat  removal  by  the  fuel  and  oxidant  flows  was 
not  provided.  Here,  convection  coefficients  were  calculated  from  an 
assumed  Nusselt  number  of  8.23  corresponding  to  fully  developed 
laminar  flow  through  an  infinitely  wide  rectangular  cross-section 
with  uniform  heat  flux  [66].  Second,  the  stated  contact  and  sepa¬ 
rator  plate  areal  resistance  was  0.1  £2  cm-2,  which  appears  to  be 
a  typographic  error  since  it  yields  an  unrealistic  90  V  ohmic  loss 
for  the  900  cm2  cell.  Rather  we  assumed  a  volume  resistivity  of 
0.1  £2  cm-3  which  yields  a  more  reasonable  0.18  V  ohmic  loss  for 
the  0.002  m  thick  interconnect  plate.  Third,  the  stated  limiting  cur¬ 
rent  was  4000  A  m-2,  but  the  I-V  curve  is  nearly  linear  and  extends 
beyond  current  densities  of  6666  Am-2.  It  was  therefore  assumed 
that  no  concentration  losses  were  present.  The  remaining  model 
parameters  were  used  as  stated  in  the  references.  The  thermal  con¬ 
ductivity  of  the  cell  and  interconnect  were  set  to  zero  since  the 
literature  cases  did  not  include  the  effects  of  in-plane  conductive 
heat  transfer. 

5.2 A.  Cases  5  and  6:  electrolyte-supported  cell 

These  models  correspond  to  the  solutions  for  an  electrolyte- 
supported  cell  solved  with  (Case  5)  and  without  (Case  6)  the 
effect  of  thermal  radiation  between  the  cell  and  interconnect.  The 
electrolyte-supported  model  consisted  of  180  p,m  electrolyte  with 
50  |jim  anode  and  cathode  electrodes.  The  cases  for  comparison 
were  performed  for  a  current  of  300  A  (3333  A  m-2 ),  fuel  utilization 
of  ~38%,  and  air  utilization  of  ~1 0%.  As  a  preliminary  test,  the  model 
was  evaluated  with  the  low  air  flow  rate  (1.09e-3  kgs-1 )  for  cur¬ 
rent  densities  of  575-6666  A  m-2  to  validate  the  implementation  of 
the  electrochemical  model.  The  results  for  the  predicted  stack  I-V 
curve  and  the  reported  reference  values  are  shown  in  Fig.  10  for 
comparison.  The  predicted  values  are  slightly  lower  than  the  refer¬ 
ence  values  with  the  largest  difference  of  -3.7%  at  the  maximum 
current  density.  As  this  depends  on  the  temperature  distributions 
through  all  of  the  cells,  the  difference  is  attributed  to  prediction 


of  a  lower  stack  temperature  distribution  as  discussed  below.  The 
electrochemical  model  is  assumed  to  correctly  capture  the  I-V  per¬ 
formance  and  therefore  validates  the  assumption  on  the  reported 
volumetric  resistivity  of  the  interconnect. 

5.2.2.  Cases  7  and  8:  anode-supported  cell 

These  models  correspond  to  the  solutions  for  an  anode- 
supported  cell  solved  with  (Case  7)  and  without  (Case  8)  the 
effect  of  thermal  radiation  between  the  cell  and  interconnect.  The 
anode-supported  cell  has  1000  fxm  anode  with  10  p,m  electrolyte 
and  25  p,m  cathode.  The  ohmic  heating  will  be  much  less  in  this 
cell  as  the  anode  resistivity  is  about  1 0,000 x  less  at  1000°C.  The 
cases  for  comparison  were  performed  for  a  current  level  of  600  A 
(6666  Am-2),  fuel  utilization  of  ~76%,  and  air  utilization  of  ~20%. 

5.2.3.  Multi-cell  stack  example  comparisons 

The  temperature  distributions  for  the  multi-cell  example  cases 
are  shown  in  Figs.  1 1  and  1 2  with  the  reported  reference  results.  The 
cell  temperature  profile  along  the  flow  direction  is  plotted  for  the 
top  and  bottom  cells  of  the  stack.  Since  the  fuel  and  oxidant  gases 
are  heated  as  they  flow  across  the  cell  due  to  the  electrochemical 
activity  in  the  cell,  the  cell  is  coolest  near  the  inlet  side  and  hottest 
at  the  outlet  side  as  expected  for  a  co-flow  configuration.  Due  to 
the  different  convective  heat  removal  capacity  between  the  fuel 
and  oxidant  streams,  an  asymmetric  temperature  profile  develops 
between  the  cells  in  the  stack,  particularly  with  an  adiabatic  bound¬ 
ary  condition.  The  low  utilization  oxidant  removes  more  heat,  so 
a  cell  at  the  periphery  of  the  stack  not  sharing  any  air  flow  chan¬ 
nel  with  adjacent  cell  will  be  cooler.  This  is  observed  in  the  results 
for  both  the  anode-supported  and  the  electrolyte-supported  where 


Fig.  10.  Comparison  of  reference  and  predicted  current-voltage  relationship  for  the 
five-cell  electrolyte-supported  stack  without  radiation. 
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Fig.  11.  Comparison  of  the  reference  and  predicted  cell  temperature  distributions  for  the  top  and  bottom  cells  of  the  five-cell  electrolyte-supported  stack  (a)  without  and 
(b)  with  radiation  heat  transfer  between  cells. 


Fig.  12.  Comparison  of  the  reference  and  predicted  cell  temperature  distributions 
with  radiation  heat  transfer  between  cells. 


the  top  and  bottom  cells  of  the  five-cell  anode-supported  stack  (a)  without  and  (b) 


the  top  cell  is  always  hotter.  This  effect  is  amplified  when  radiation 
between  cells  is  neglected.  Without  radiation,  convection  between 
the  solid  and  fluid  domains  is  the  only  heat  transfer  mechanism 
from  cell  to  cell.  With  radiation,  dual  heat  transfer  mechanisms 
reduce  the  temperature  difference  between  cells  and  drive  the 
stack  to  a  more  isothermal  condition.  This  was  demonstrated  for 
both  cell  support  types  and  indicates  the  importance  of  radiation 
for  reducing  thermal  gradients. 

Comparison  of  the  results  shows  that  the  trends  are  well  cap¬ 
tured  but  some  important  differences  between  the  temperature 
fields  remain.  Most  noticeably,  the  prediction  for  the  electrolyte- 
supported  cell  without  radiation  is  low  overall.  Other  differences 
include  a  generally  higher  prediction  for  the  anode-supported  cell 
with  radiation  and  a  different  curve  shape  for  the  anode-supported 
cell  without  radiation.  The  differences  can  be  attributed  primarily 
to  the  assumptions  on  the  convective  heat  transfer.  First,  because 
the  empirical  Nusselt  number  correlation  is  not  known,  the  use 
of  the  relation  for  a  parallel  plate  enclosure  will  introduce  some 
difference.  Second,  the  model  in  its  current  state  assumes  only  a 
single  value  for  convection  coefficient  along  the  flow  domain  and 
neglects  any  temperature  effects  for  both  the  oxidant  and  fuel. 
Furthermore,  the  convective  heat  transfer  coefficient  for  the  fuel 
may  also  depend  on  the  composition,  as  the  thermal  conductiv¬ 
ity  of  hydrogen  is  4x  higher  than  steam  at  1200K  (528e-3  and 
1 25e-3  W  K-1  m-1 ,  respectively).  The  convection  coefficient  of  the 
mixture  will  vary  from  the  fuel  rich  inlet  to  the  fuel  depleted  out¬ 
let.  The  model  here  used  a  single  value  based  on  the  average  cell 
temperature  and  average  composition,  so  this  is  assumed  to  be 


the  primary  assumption  responsible  for  the  difference  between 
the  results  predictions.  Because  lower  temperature  yields  higher 
Ohmic  loss,  a  slightly  lower  I-V  response  is  predicted  as  shown 
in  Fig.  10.  Overall,  reasonably  good  model  comparisons  have  been 
achieved. 

6.  Simulation  results  for  multi-cell  stacks 

The  SOFC-MP  modeling  tool  was  next  used  to  evaluate  a  large 
multi-cell  stack  characteristic  of  the  modern  stacks  currently  under 
development  in  the  SECA  program.  This  model  served  as  a  baseline 
for  comparison  with  other  models  that  demonstrate  the  features 
of  this  tool.  Specifically,  the  tool  was  used  to  evaluate  the  cooling 
benefits  of  a  reforming  fuel.  The  model  also  has  the  capability  for 
simulating  the  cell-to-cell  variations  with  different  stack  param¬ 
eters.  This  feature  is  useful  for  identifying  the  sensitivity  of  stack 
performance  to  various  off-normal  conditions,  so  that  researchers 
can  maximize  their  understanding  of  stack  performance  in  actual 
tests  and  experiments.  This  capability  for  cell-to-cell  variation 
which  is  absent  in  most  0D  and  ID  models  and  proven  computa¬ 
tionally  expensive  for  3D  models  can  be  accomplished  fairly  easily 
in  the  SOFC-MP  2D  model.  The  model  is  currently  capable  of  simu¬ 
lating  the  following  cell-to-cell  variations,  e.g.: 

•  Thick  plates  in  middle  of  stack  (e.g.,  for  thermocouple  measure¬ 
ments,  load  frame). 

•  Different  flow  rates  in  cells  (e.g.,  blockage,  leak,  by-pass). 
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Fig.  13.  High  performance  I-V  response  for  the  baseline  96-cell  stack  model  with 
both  hydrogen  and  reforming  fuels. 

•  Different  I-V  performance  of  cells. 

•  Short  current  in  cells. 

•  Partial  contact  in  cells. 

As  examples  in  this  paper,  the  influence  of  an  instrumented  plate 
introduced  to  the  stack  for  temperature  measurements  and  a  gas 
flow  restriction  is  assessed.  In  this  section,  the  stack-wise  metric  for 
simulation  of  the  96-cell  stack  is  obtained  for  comparison  to  several 
off-normal  conditions  to  demonstrate  the  capabilities  of  SOFC-MP 
2D  model. 

6. 1.  96-Cell  baseline  stack 

SOFC  manufactures  have  made  great  technical  strides  in  the  pur¬ 
suit  of  larger,  more  powerful  stacks.  Developers  within  the  SECA 
program  are  currently  testing  stack  towers  with  normal  power 
outputs  of  25  kW  or  larger  [67].  For  these  stack  towers  consisting 
of  nearly  200  planar  cells,  variation  of  temperature  and  electrical 
performance  is  expected  along  the  stack  height.  These  variations 
must  be  understood  to  determine  appropriate  operating  conditions 
that  prevent  high  temperatures,  high  current  densities,  high  local 
fuel  utilizations,  or  other  behaviors  that  can  lead  to  performance 
degradation.  In  this  section,  the  SOFC-MP  2D  model  will  be  used  to 
evaluate  flow-thermal-electrochemical  performance  of  a  tall  stack 
characteristic  of  modern  designs  while  still  providing  resolution  of 
the  field  variation  across  cells  and  along  the  stack  height.  The  simu¬ 
lations  are  based  on  a  hypothetical  model  due  to  the  lack  of  detailed 
published  data  on  tall  SOFC  stacks. 

The  baseline  stack  consisted  of  a  large  625  cm2  active  area 
cell  (25  cm  x  25  cm)  operated  at  an  average  current  density  of 
0.5  A  cm-2  in  a  96-cell  stack  to  provide  a  total  power  output  p 

of  25.5  kW.  The  temperature  of  the  inlet  gases  was  assumed  to 
increase  along  the  stack  height  from  700  °C  for  the  bottom  cell 
to  730  °C  at  the  top  cell  of  the  stack.  The  stack  was  assumed  to  2 

be  in  a  750  °C  furnace  with  external  stack  heat  loss  via  radiation  ^ 

(emissivity  of  0.7)  and  convection  (free  convection  coefficients  of  £ 

0.3-1 .5  W m-2 1<).  The  electrochemistry  parameters  were  selected  ® 

to  give  a  very  high  performing  cell  characteristic  of  the  state-of- 
art  for  peak  power  conditions  of  SECA  SOFCs  [38,39,58],  and  the 
I-V  curve  is  shown  in  Fig.  13.  The  fuel  was  wet  hydrogen  with  3% 
water,  and  the  stack  was  run  at  65%  fuel  utilization  and  15%  air 
utilization. 

The  predicted  performance  was  an  average  cell  voltage  of 
0.852  V  and  total  power  of  25.5  kW.  Thermal  performance  indicated 
a  mean  cell  temperature  of  793  °C  with  minimum  and  maximum 


Fig.  14.  Cell  voltage  profile  for  the  96-cell  stack  operating  with  wet  hydrogen  fuel. 


Cell  position  /  m 


Fig.  15.  Current  density  distribution  along  various  cells  of  the  96-cell  stack  operat¬ 
ing  with  wet  hydrogen  fuel. 

values  of  709  °C  and  848  °C,  respectively.  Results  for  the  voltage 
profile  along  the  stack  height,  current  density  across  various  cells, 
cell  peak  temperatures,  cell  temperature  across  various  cells,  and 
cell  temperature  difference  is  shown  in  Figs.  14-18,  respectively. 
The  profile  in  Fig.  14  showed  that  the  voltage  generally  increased 
along  the  stack  height  with  a  total  variation  of  0.015  V,  while  the 
top  and  bottom  cells  showed  a  local  voltage  decrease.  This  per¬ 
formance  is  due  to  the  improved  electrochemical  performance  for 


Fig.  16.  Profile  of  cell  minimum,  maximum,  and  average  cell  temperatures  for  the 
96-cell  stack  operating  with  wet  hydrogen  fuel. 
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Fig.  17.  Temperature  distribution  along  various  cells  for  the  96-cell  stack  operating 
with  wet  hydrogen  fuel. 

regions  with  higher  temperature.  The  temperature  profile  (Fig.  16) 
shows  that  the  cells  are  generally  hotter  near  the  top  of  the  stack 
due  to  the  inlet  gas  flow  temperatures,  but  heat  transfer  with  the 
furnace  affects  the  temperatures  for  the  top  and  bottom  cells  by 
driving  them  toward  750  °C.  The  current  density  for  the  top,  bot¬ 
tom,  and  middle  cells  are  shown  in  Fig.  15.  Cells  higher  in  the  stack 
tend  to  consolidate  the  current  near  the  fuel  inlet  with  the  hotter 
inflow  gases.  Flowever,  the  minimum  temperature  of  the  top  and 
bottom  cells  (which  occurs  at  the  fuel  inlet  as  shown  in  Fig.  17) 
increases  due  to  furnace  heating  so  that  these  end  cells  show  the 
most  current  consolidation  at  the  leading  edge  (Fig.  15).  This  will 
be  important  to  stack  designs  since  high  current  densities  are  likely 
to  be  more  conducive  to  various  degradation  processes.  Also,  the 
highest  cell  temperature  difference  of  about  139  °C  occurred  in  cell 
#15  (Fig.  18).  The  cell  temperature  difference  is  also  an  important 
metric  for  design  since  greater  thermal  strain  variation  can  lead  to 
unwanted  cell  warpage  and  structural  stresses. 

6.1.1.  Effect  of  reforming 

The  use  of  on-cell  reforming  of  methane  is  advantageous  for 
SOFCs  because  the  endothermic  reaction  removes  heat  and  may 
be  beneficial  in  thermal  management  for  cells  with  large  area.  The 
baseline  model  was  run  at  the  same  current  density  with  a  50% 
on-cell  reforming  fuel  consisting  of  a  molar  composition  of  32.4% 
H2,  33.3%  H20,  4.9%  CO,  6.1%  C02,  11.0%  CH4  and  12.4%  N2.  The 
reforming  rate  following  the  expression  of  Achenbach  was  used, 


Fig.  18.  Cell  temperature  difference  profile  for  the  96-cell  stack  operating  with  wet 
hydrogen  fuel. 


Fig.  19.  Voltage  profile  for  the  96-cell  with  reforming  fuel. 


Fig.  20.  Profile  of  the  minimum,  maximum,  and  average  cell  temperatures  for  the 
96-cell  with  reforming  fuel. 

and  the  CFI4  was  fully  reformed  for  this  case.  The  total  power  out¬ 
put  decreased  7.5-23.6  kW  due  to  the  lower  Nernst  potential  which 
resulted  in  an  average  cell  voltage  of  0.785  V.  The  voltage  profile  in 
Fig.  19  shows  that  the  lowest  voltage  was  observed  near  cell  #15. 
This  was  again  explained  by  the  temperature  profile  and  the  periph¬ 
ery  boundary  condition  effects  which  lead  to  the  lowest  average  cell 
temperature  to  occur  in  this  part  of  the  stack  (Fig.  20).  The  current 
density  distribution  (Fig.  21 )  shows  the  highest  values  near  the  cell 
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Fig.  21.  Current  density  distribution  for  various  cells  of  the  96-cell  stack  with 
reforming  fuel. 
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Fig.  22.  Demonstration  of  the  effect  of  using  thicker  interconnect  plates  for  cells  #30  and  #60  on  the  temperature  (a)  peak  values  and  (b)  distribution. 


outlet  due  to  the  endothermic  reaction  near  the  fuel  inlet  side.  The 
peak  current  density  in  the  stack  was  increased  6-0.69  A  cm-2  com¬ 
pared  to  the  peak  value  of  0.65  A  cm-2  in  the  H2  stack.  The  benefits 
of  the  reforming  fuel  were  observed  when  the  average  cell  tempera¬ 
ture  was  reduced  from  793  °C  to  743  °C,  the  maximum  temperature 
reduced  31  °C  from  848  °C  to  817°C,  and  the  maximum  cell  tem¬ 
perature  difference  reduced  20  °C  from  139  °C  to  119°C.  Models 
such  as  these  can  help  the  stack  designer  with  the  judicious  use  of 
reforming  to  manage  trade-offs  between  high  temperatures  that 
can  lead  to  material  degradation,  high  temperature  gradients  that 
can  create  thermal-mechanical  stresses,  and  high  current  densities 
that  can  lead  to  degradation. 

6.1.2.  Effect  of  measurement  plates 

Fuel  cell  developers  have  sought  ways  to  obtain  real  time  mea¬ 
surement  of  the  stack  thermal  state  during  testing  to  assess  the 
design  margin  between  maximum  operating  temperature  and  tem¬ 
perature  limits  for  material  stability.  This  is  critical  in  achieving  an 
acceptable  long-term  power  degradation  rate  by  preventing  var¬ 
ious  degradation  mechanisms  that  are  thermally  activated  (e.g., 
scale  growth,  oxidation,  interfacial  reactions,  species  volatilization, 
etc.).  Evaluation  of  temperature  within  the  stack  is  difficult  because 
precise  sensor  placement  will  likely  require  an  undesirable  pene¬ 
tration  through  the  hermetic  seals  in  the  stack,  and  multiple  sensor 
placements  within  manifolds  and  cells  may  adversely  affect  the  gas 
flow  characteristics,  and  any  instrumentation  must  ensure  it  does 
not  provide  a  short  circuit  electrical  path.  One  effective  alterna¬ 
tive  is  to  use  an  instrumented  plate  which  consists  of  numerous 
embedded  thermocouples  in  a  structure  patterned  as  the  cell  man¬ 
ifold  design  which  can  then  be  introduced  and  sealed  between  any 
two  cells  in  a  tall  stack.  The  numerous  sensors  can  then  provide 
enough  resolution  to  characterize  the  actual  temperature  profile 
across  the  cell. 

However,  this  approach  could  potentially  suffer  from  the  so- 
called  “observer  effect”  where  the  measurement  method  strongly 
influences  the  data  being  measured.  The  thicker  instrumented  plate 
will  likely  provide  an  enhanced  in-plane  thermal  conduction  com¬ 
pared  to  the  thin  interconnect  which  will  affect  the  measured 
temperature  field.  In  this  example,  the  plate  thickness  of  cells  #32 
and  cell  #64  is  increased  from  the  nominal  value  of  0.5-12  mm, 
simulating  the  introduction  of  two  relatively  thick  measurement 
plates.  Fig.  22(a)  shows  the  maximum,  minimum,  and  average  tem¬ 
perature  distribution  of  all  cells  in  the  stack,  while  Fig.  22(b)  plots 
the  temperature  distribution  in  the  flow  direction  for  measured  cell 
#30  and  its  neighbor  cell  #32.  Comparison  of  Fig.  22  to  the  nominal 
case  in  Fig.  1 6  shows  that  while  the  average  temperature  of  the  cells 


is  nearly  the  same  in  the  neighborhood  of  the  measurement  plates, 
the  minimum  and  maximum  temperatures  of  these  cells  are  greatly 
affected.  For  the  thicker  plates,  the  improved  thermal  conductiv¬ 
ity  path  spreads  heat  better  to  smooth  temperature  gradients  in 
the  cell.  The  minimum  temperature  at  the  inlet  is  higher  while  the 
maximum  temperature  is  lower.  This  results  in  the  designer  under¬ 
estimating  the  temperature  difference  by  18%  and  the  maximum 
temperature  by  13  °C.  This  underestimation  is  highly  undesirable 
since  it  may  mislead  the  designer  to  think  that  the  stack  is  operating 
within  the  desired  temperature  range  while,  in  fact,  it  is  not.  Similar 
thermal  effects  are  expected  from  the  use  of  large  manifold  or  load 
frame  structures  in  tall  stack  towers.  Modeling  tools  such  as  this 
can  help  the  designer  to  understand  the  impact  of  these  structures 
on  the  thermal  profile  in  the  stack. 

6. 1.3.  Effect  of  flow  maldistribution 

The  cells  in  the  SOFC  system  are  generally  designed  to  run  at  an 
optimal  condition  with  uniform  flow  rates  in  each  cell  to  achieve 
the  desired  fuel  and  air  utilization.  However,  certain  design  and 
assembly  factors  can  adversely  affect  gas  flows  (e.g.,  leak  through 
poor  seals,  flow  restriction/blockage  due  to  sealant  overflow,  by¬ 
pass  of  fuel  around  the  cell)  and  downgrade  the  cell  performance. 
The  SOFC-MP  2D  model  can  be  used  to  evaluate  how  localized 
downgraded  cells  will  impact  the  overall  stack  performance.  As  an 
example  case,  cells  #32  and  #64  of  the  96-cell  stack  are  assumed 
to  have  flow  restrictions,  i.e.,  cell  #32  has  a  25%  reduction  in  fuel 
flow  rate  and  cell  #64  has  a  50%  reduction  in  oxidant  flow  rate. 

The  results  show  that  the  fuel  flow  restriction  has  the  most  sig¬ 
nificant  effect  on  the  performance.  The  reduced  fuel  flow  increases 
the  fuel  utilization  on  cell  #32  from  65%  in  the  baseline  case  to  87%. 
This  reduces  the  H2  concentration  over  the  cell  which  reduces  the 
Nernst  potential,  makes  cell  #32  less  efficient,  and  drops  its  voltage 
by  5%  from  the  baseline  0.85-0.81  V.  The  current  then  concentrates 
near  the  fuel-rich  inlet  with  a  21%  increase  in  peak  current  den¬ 
sity  from  0.61  A  cm-2  in  the  baseline  case  to  0.74  A  cm-2  (Fig.  23a). 
The  higher  heat  load  then  increases  the  average  local  cell  temper¬ 
ature  by  9°C  (+1.1%)  but  the  peak  cell  temperature  only  increases 
by  5°C  (+0.6%)  as  shown  in  Fig.  23b.  The  flow  restriction  on  the 
air  had  a  much  less  impact  since  the  air  utilization  is  much  lower. 
The  voltage  decreased  slightly  by  0.5  V  from  the  adjacent  cells,  but 
the  effect  of  air  restriction  on  the  cell  temperature  was  almost 
negligible. 

6.1.4.  Effect  of  electrochemical  performance 

The  variation  of  cell  performance  depends  not  only  on  the  oper¬ 
ating  conditions  but  also  on  the  quality  control  of  the  cell  fabrication 


3220 


K.  Lai  et  al.  /  Journal  of  Power  Sources  196(2011)  3204-3222 


Fig.  23.  Contours  of  the  (a)  minimum,  maximum,  and  mean  cell  temperatures  and 
flow  rates. 
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current  density  distribution  for  the  96-cell  stack  with  reduced  fuel  (#32)  and  air  (#64) 
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Fig.  24.  Voltage  and  temperature  distribution  for  96-cell  stack  with  four  poor  performing  cells:  #32,  #33. 


and  assembly  processes.  Cells  with  inherently  different  electro¬ 
chemical  properties  may  be  present  in  one  single  assembled  stack. 
Alternatively,  undesirable  degradation  processes  may  occur  within 
a  tested  stack  that  cause  local  voltage  drops.  The  modeling  tool  can 
simulate  the  different  electrochemical  performance  of  individual 
cells  to  simulate  the  effect  of  underperforming  cells  on  the  stack 
response.  This  is  important  as  it  provides  the  designer  with  the 
maximum  amount  of  information  in  a  stack  experiment  where  volt¬ 
age  degradation  has  occurred.  For  this  example,  four  cells  #32-35 
were  assumed  to  have  added  resistance  resulting  in  an  additional 
0.1  V  voltage  loss. 

Fig.  24(a)  shows  the  stack-wise  voltage  distribution  and  (b) 
shows  the  temperature  distribution  changes  caused  by  the  down¬ 
graded  cells.  The  stack  power  output  correspondingly  decreased 
about  0.35%,  but  the  primary  influence  on  the  stack  was  local  tem¬ 
perature  increase  in  the  cells.  The  added  ohmic  heat  increased  the 
local  average  cell  temperature  by  ~6  °C,  the  peak  cell  temperature 
by  ~7°C,  and  the  cell  temperature  difference  by  ~6°C.  Depend¬ 
ing  on  the  severity  of  voltage  degradation  and  the  number  of  cells 
affected,  the  impact  on  the  stack  temperatures  may  be  important 
for  certain  components’  temperature  limits. 

7.  Model  performance  and  scalability 

The  total  number  of  control  volumes  in  the  system  is 
(Nceiis  x  4  +  1)NX  4 NceUsNx,  which  is  proportional  to  the  number 
of  cells  Ncells  and  the  number  of  increments  in  the  axial  direction 
Nx.  A  large  portion  of  the  computational  time  is  spent  on  solving 
the  matrix  for  temperature.  The  computational  time  spent  on  the 


matrix  solver  chosen  for  this  sparsely  populated  thermal  matrix  is 
slightly  more  than  linear  to  the  total  unknowns.  The  total  itera¬ 
tions  from  the  two  loop  iterations  (temperature  field  and  current) 
depends  more  on  the  multi-physics  characteristics  of  the  model 
rather  than  the  total  number  of  control  volumes.  Therefore,  the 
total  computation  time  is  proportional  to  the  number  of  cells  and 
the  number  of  mesh  increments.  Naturally,  the  computational  time 
also  depends  on  the  convergence  tolerance  set  for  the  temperature 
field  and  the  current. 

It  is  shown  that  for  all  benchmark  cases,  mesh  size  convergence 
is  usually  achieved  around  Nx  =  200,  regardless  of  the  number  of 
cells  or  the  flow  direction.  As  a  result,  the  SOFC-MP  model  can  eas¬ 
ily  manage  the  simulation  of  a  large  system,  e.g.,  a  system  with 
100  cells  and  500  increments.  For  co-flow,  a  system  with  100  cells 
and  500  increments  requires  about  8  min  on  a  2  GHz  Windows 
XP  server  without  the  restart  capability  using  the  previous  saved 
data.  Because  the  air  flows  in  the  opposite  direction,  the  thermal 
properties  of  the  air  will  not  be  known  from  the  previous  control 
volume;  they  will  be  assigned  values  from  previous  iteration  step 
instead.  As  a  result,  it  takes  more  iteration  steps  to  achieve  a  stable 
solution,  and  the  computational  time  on  a  counter-flow  system  is 
usually  higher,  e.g.,  about  20  min  for  the  same  configuration,  but 
still  manageable. 

The  computational  effort  for  SOFC-MP  simulation  is  much 
greater  if  the  average  current  density  solution  mode  is  specified 
instead  because  multiple  simulations  for  the  trial  average  voltage 
will  be  run.  A  typical  simulation  for  target  current  density  runs 
about  10  times  the  equivalent  simulation  on  voltage,  but  is  still 
manageable. 
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8.  Conclusions 

A  quasi-2D  multi-physics  finite  volume  SOFC  model  was 
successfully  developed  for  simulating  steady-state  co-flow  and 
counter-flow  multi-cell  SOFC  stack  systems.  Species  conservation, 
energy  conservation  equations,  and  electrochemical  I-V  models 
were  solved  in  the  model  to  provide  the  overall  system  perfor¬ 
mance  metrics  along  with  the  distributions  of  temperature,  species 
concentrations,  current  density,  and  voltage  through  the  stack.  The 
predicted  results  from  the  model  compared  very  well  with  sin¬ 
gle  cell  benchmark  cases  using  different  fuel  compositions  and 
flow  geometries.  The  results  for  five-cell  co-flow  and  counter¬ 
flow  stacks  with  and  without  internal  radiation  also  compared 
reasonably  well  with  published  simulation  results.  The  authors  rec¬ 
ommend  that  a  comprehensive  suite  of  well  defined  benchmarks 
should  be  established  by  the  fuel  cell  community  to  support  SOFC 
model  validation  activities  and  include  potential  variations  for  cell 
size,  stack  size,  external  heat  loss,  on-cell  reforming,  and  manifold 
structures. 

Simulation  results  have  been  presented  for  a  25  kW  multi-cell 
stack  to  illustrate  the  model  capability  on  handling  cell-to-cell 
variations  and  user  defined  electrochemistry.  On-cell  reforming 
was  shown  to  beneficially  reduce  the  stack  temperature  and  tem¬ 
perature  gradient  at  the  expense  of  a  small  power  loss.  The 
presence  of  thick  plates  in  the  stack  for  thermocouple  placement 
was  shown  to  dramatically  affect  the  stack’s  measured  peak  tem¬ 
perature  and  temperature  difference,  indicating  that  designers 
must  included  these  features  to  accurately  characterizes  multi¬ 
cell  stack  experiments.  The  influence  of  reduced  fuel/oxidant 
flow  and  cell  electrochemical  performance  on  the  maximum 
current  density  and  cell  temperatures  was  also  demonstrated. 
Future  study  will  include  a  more  comprehensive  parametric 
study  on  large  cell  stacks  using  the  cell-to-cell  variation  fea¬ 
ture  to  characterize  stack  performance  for  realistic  performance 
variations. 
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